An upper-crust lid over the Long Valley magma chamber

Geophysical characterization of calderas is fundamental in assessing their potential for future catastrophic volcanic eruptions. The mechanism behind the unrest of Long Valley Caldera in California remains highly debated, with recent periods of uplift and seismicity driven either by the release of aqueous fluids from the magma chamber or by the intrusion of magma into the upper crust. We use distributed acoustic sensing data recorded along a 100-kilometer fiber-optic cable traversing the caldera to image its subsurface structure. Our images highlight a definite separation between the shallow hydrothermal system and the large magma chamber located at ~12-kilometer depth. The combination of the geological evidence with our results shows how fluids exsolved through second boiling provide the source of the observed uplift and seismicity.


Supplementary Text
Checkerboard test and model resolvability We conduct conventional seismic tomography checkerboard tests in which oscillatory anomalies of -5% and 5% variations in the P-(VP) and S-wave (VS) speeds are introduced within the initial models (40,68).Based on the estimated picking errors, we introduce zero-mean Gaussian noise to the model traveltimes whose standard deviations are 50 ms and 100 ms to the P-and S-wave picks, respectively.We also randomly perturb the event locations by ±500 m in each direction.The panels in Fig. S3 show the true and retrieved perturbations at various depth levels for the inverted VS.To indicate the portions of the models that are resolvable by the tomography approach, we employ the resolvability index (5,69,70).This index ranges from 0 to 1: 0.5 represents portions of the model for which 0% of the perturbation is retrieved (i.e., zero sensitivity), a value of 1 is a perfect reconstruction, and 0 defines portions for which -100% of the perturbation is inverted.The shaded areas in these panels indicate a resolvability index smaller than 0.6, which is considered a reasonable lower bound for resolvable areas (5,70).Similar results are obtained for the VP perturbations.The best resolution (i.e., no smearing effect) is obtained in the proximity of the DAS channels for most of the depth slices.For the shallow slices (-2 to 2 depths), almost no smearing is observed in the proximity to the sensed fibers.Deeper slices present a broader area of resolvability but at the cost of smearing the perturbations.Below 15 km depth, a significantly reduced portion of the subsurface is resolvable due to the limited number of rays reaching that portion of the model.P-wave velocity anomalies, residual histograms, and velocity-model validation Fig. S4 displays the same slices of Fig. 2 in the main text but for the VP anomalies from a onedimensional Walker Lane crust velocity model shown in Fig. S5.This one-dimensional model is obtained by averaging the initial model along the latitude and longitude axes.Compared to the initial model shown in the panels on the left column, the inverted velocity anomalies (right column) present the same clear separation between the magmatic chamber centered at approximately 12.5 km depth and the shallow structures above as in the S-wave anomalies.Additionally, similarly to the VS model, velocity reductions within the caldera, along the Mono-Inyo craters, and underneath the Mono lake are obtained by the tomographic workflow.Similar features are obtained from the tomographic workflow when using this 1D model as an initial guess.However, these anomalies are better defined by starting the tomographic process from the surface-wave inverted velocities.In addition, the 3D surface-wavederived model as an initial guess provides lower final traveltime residuals compared to the 1D model when used to start the tomography workflow.
The top panels in Fig. S6 show the P-and S-wave absolute traveltime residual histograms obtained using the initial velocity models (Fig. S6 A and B, respectively).Their respective residual means are -0.17s and 0.11 s, while their standard deviations are 0.65 s and 0.78 s.The relocation and tomography workflow produces velocity models whose traveltime residuals are Gaussian distributed with means of -0.01 and 0.01 and standard deviations of 0.4 s and 0.47 s for P-and S-wave traveltimes, respectively (Fig. S6 C and D).Tighter distributions could be achieved by relaxing the smoothness constraints defined by the Gaussian filter employed during the inversion process.
However, smaller-scale velocity anomalies are not resolvable by the event-channel geometry (based on checkerboard test analyses) and thus we avoid introducing them during the inversion process.
Finally, we validate our inverted velocities by modeling the arrival times for a relocated event that was not included during the tomographic steps.The event ID from the NCEDC DD catalog is 73485976 and its magnitude and relocated depth are 2.8 and 2.327 km, respectively.The maximum and minimum distances from the DAS channels are 33.5 km and 54.5 km (Fig. S7).Fig. S7 displays the recorded DAS data in which we overlay the P-(red lines) and S-wave (blue lines) traveltimes predicted from the initial (dashed lines) and inverted (solid lines) velocity models.The final velocities predict traveltimes closely following the observed arrival onsets.On the other hand, the initial models underestimate the observed onsets, highlighting the presence of low-velocity anomalies (e.g., Mono-Inyo crater basin) necessary to obtain correct traveltime predictions.
Our melt fraction estimations are based on a linearized VS/melt-fraction derivative (δVS/δMF) of -23 m/s/MF derived by averaging the Voigt and Reuss VS/melt-fraction trends for a 4% H2O-wet rhyolite at 310 MPa and 750 Cº (27).In our models, the -15% and -20% S-wave reductions within the magma chamber correspond to wave speeds of 3.0 km/s and 2.86 km/s with VP/VS of 1.83 and 1.86.

Hypothesis-driven resolution tests
To test the resolution and bias of our workflow in imaging known seismic anomalies, we invert synthetic traveltime datasets where different type of velocity reductions.As in other tomographic studies (5,32), when displaying the results of our synthetic tests, we show the input and inverted anomalies with the background model removed to better highlight the resolved portion of the introduced perturbation.
In our first resolvability test, we introduce a low VP and VS anomaly underneath Mono Lake with a high VP/VS to simulate the presence of a large water saturated materials (Fig. S8A-D).We apply the same procedure as described for the checkerboard test to invert the synthetic traveltimes.Our tomographic workflow can correctly retrieve the shape and the overall velocity decrease as well as VP/VS ratio with a minor underestimation of both properties (Figs.S7E-H).
In our second set of simulated experiments (Figs.S8-S10), we incorporate a sequence of lowshear-wave velocity reductions that progressively diminish in size.The goal is to evaluate the resolution limits of our method in imaging upper-crust magma reservoirs, which are not present within our results obtained using our DAS dataset.In each test we introduce a cylinder-shaped anomaly, a large perturbation of 10 km radius and 4 km thickness (Fig. S9), a middle-size reservoir of 5 km radius and 4 km thickness (Fig. S10), and a small 2 km radius and 2 km thickness anomaly (Fig. S11).In all three cases the shape of the anomalies is well-imaged by our tomographic approach with an underestimation of the velocity decrease within the center of the anomaly; especially, for the smallest anomaly of 2 km radius.If the large and middle-size reservoirs were present within our results, we would have clearly detected and interpreted them as potential upper-crust magma reservoirs.On the contrary, the smallest anomaly is close to the resolution limits of the method given the employed DAS geometry and earthquake locations.Thus, any small upper-crust reservoirs whose core has a size less than 2 km in diameter may be challenging to detect due to the underestimation of the velocity reduction.

Fig. S1 .
Fig. S1.Local and regional events used in this study.Conventional stations are depicted by the blue triangles while the green line indicates the locations of the DAS channels.The earthquakes are indicated by the red dots whose sizes are proportional to their magnitude.

Fig. S2 .
Fig. S2.Assessing picking accuracy using event cross-correlations.(A) Crosscorrelation example between two events recorded by the South DAS array.The cyan line represents the cross-correlation-based shift retrieved by the multi-channel crosscorrelation algorithm.(B and C) Differential traveltime histograms for the P-and S-wave cross-correlation windows, respectively.

Fig. S3 .
Fig. S3.Checkerboard test results for the VS model.The size of each Gaussian anomaly is 10 km in each direction.

Fig. S4 .
Fig. S4.The Long Valley P-wave anomalies.The panels on the left column display the initial model structures, while the panels on the right depict the inverted P-wave velocity anomalies.All perturbations are with respect to a one-dimensional Walker Lane crust profile (obtained by averaging the initial model along latitude and longitude).(A to B) Depth slices at -1.0 km elevation.The caldera and lakes' extents are shown by the black dashed lines.(C to H) Model profiles indicated in panel A. The white (panels A and B) and black (panels C to H) dashed lines indicate the -15% and -12% P-wave velocity contours.The white solid lines separating the shaded areas denote the resolvable model portions.

Fig. S5 .
Fig. S5.Reference Walker Lane velocity profiles.P-and S-wave speed profiles employed to compute all the velocity perturbations shown in Figs. 2 and S4.

Fig. S6 .
Fig. S6.Initial and final traveltime residual histograms.(A and B) Absolute traveltime residual histograms obtained using the initial P-and S-wave speed models.(C and D) Same as the top panels but obtained with traveltimes predicted using the inverted models.

Fig. S7 .
Fig. S7.Inverted wave-speed assessment with an independent event.(Left) Map showing the location of the DAS array channels (black solid line) and of the event (red star) used for result validation.Recorded DAS strain of the event overlaid by the P-(red lines) and S-wave (blue lines) predicted arrival times using the initial (dashed lines) and inverted (solid lines) models (Right).

Fig. S8 .
Fig. S8.Resolvability test for low velocities and high VP/VS anomaly under Mono Lake.(A-D) VS decrease and VP/VS ratio anomaly introduced within our model to test the resolvability of the tomographic workflow underneath Mono Lake.A similar synthetic VP reduction is placed under Mono Lake.(E-H) VS reduction and VP/VS ratio anomaly retrieved by our method.The arrows point at the introduced and recovered anomaly in each panel.The black dashed lines in panels B and F indicate -15% and -20% velocity reductions.

Fig. S9 .Fig. S10 .
Fig. S9.Resolvability test for an upper-crust low-velocity anomaly under the Long Valley caldera of 10 km radius.(A) Input synthetic VS anomaly mimicking an uppercrust magma reservoir.(B) Velocity anomaly retrieved by our tomographic workflow.The black dashed lines in the velocity profiles indicate -15% and -20% velocity reductions.